import os
import warnings
import subprocess
from pathlib import Path
import numpy as np
import pandas as pd
from sklearn.decomposition import PCA
import seaborn as sns
import plotly.express as px
import matplotlib.pyplot as plt
sns.set(style="whitegrid", context="notebook", palette="pastel")
warnings.simplefilter('ignore')
DATA_FOLDER = Path('../data/HW1')
SCRIPT_FOLDER = Path('../scripts')
SEED = 42
🧬 Task 1¶
Для набора из трех позиций и четырех сэмплов, сделайте PCA. Спроектируйте на главные компоненты, мотивировав чем-то выбор количества компонент.
Матрица:
$$X = \begin{pmatrix} 0 & 1 & 0 \\ 0 & 0 & 1 \\ 1 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix} $$
Решение:
- Центрируем матрицу $$ Z = X - \bar{X} $$
- Средние по столбцам: $\bar{X_1} = \bar{X_2} = \bar{X_3} = 0.25$
- Стандартизировать не будем, все и так в одном масштабе
$$ Z = \begin{pmatrix} -0.25 & 0.75 & -0.25 \\ -0.25 & -0.25 & 0.75 \\ 0.75 & -0.25 & -0.25 \\ -0.25 & -0.25 & -0.25 \\ \end{pmatrix} $$
- Строим ковариационную матрицу
Ковариация (несмещенная оценка): $$ \text{cov}(X, Y) = \frac{1}{n-1} \sum_i^n (X - \bar{X}) \cdot (Y - \bar{Y}) $$
Ковариационная матрица - матрица, у которой $(i,j)$-элемент является ковариацией признаков $(X_i, X_j)$. Когда $X_i = X_j$ ковариация - это дисперсия $X_j$ $(\mathbb{D}X_j)$
- По диагонали - дисперсии признаков
- В остальных ячейках - ковариации соответствующих признаков
- Т.к. ковариация симметричная, то и матрица ковариации симметрична
$$ \text{COV\_MATRIX} = \begin{pmatrix} \text{cov}(X_1, X_1) & \text{cov}(X_1, X_2) & \text{cov}(X_1, X_3) \\ \text{cov}(X_2, X_1) & \text{cov}(X_2, X_2) & \text{cov}(X_2, X_3) \\ \text{cov}(X_3, X_1) & \text{cov}(X_3, X_2) & \text{cov}(X_3, X_3) \\ \end{pmatrix} $$
Наша ковариационная матрица: $$ \mathbb{D}X_1 = \mathbb{D}X_2 = \mathbb{D}X_3 = \frac{(1 + 1 + 1 + 9)}{3 \cdot 4^2} = \frac{1}{4} $$
$$ \text{cov}(X_1, X_2) = \text{cov}(X_1, X_3) = \text{cov}(X_2, X_3) = \frac{-3 + 1 - 3 + 1}{3 \cdot 4^2} = -\frac{1}{12} $$
$$ \text{COV} = \begin{pmatrix} \frac{1}{4} & -\frac{1}{12} & -\frac{1}{12} \\ -\frac{1}{12} & \frac{1}{4} & -\frac{1}{12} \\ -\frac{1}{12} & -\frac{1}{12} & \frac{1}{4} \\ \end{pmatrix} $$
- Ищем собственные числа и вектора в ковариационной матрице.
Собственные значения находят через уравнение (3), после чего с известными с.з. находят собственные вектора через уравнение (2) $$ \begin{gather} A \cdot v = \lambda \cdot v \\ v \cdot (A - \lambda I) = 0 \\ det(A - \lambda I) = 0 \end{gather} $$ где $A$ - исходная матрица, $v$ - собственный вектор, $\lambda$ - собственное значение, $I$ - identity matrix
- Каждому собственному вектору соответствует свое собственное значение (пары собственные вектор-значение)
- Собственные значения отражают "важность" компоненты - какую долю дисперсии изначальных признаков она несет (больше дисперсия $\to$ больше информации содержит данный признак. Т.е. чем больше собственное значение, тем важнее компонента)
- Каждый собственный вектор показывает направление, вдоль которого наблюдается соответствующая дисперсия. Это и есть компоненты в новом пространстве признаков!
- Компоненты в новом пространстве независимы друг от друга
$$ \begin{gather} det(A - \lambda I) = 0 \\ \begin{vmatrix} \frac{1}{4} - \lambda && -\frac{1}{12} && -\frac{1}{12} \\ -\frac{1}{12} && \frac{1}{4} - \lambda && -\frac{1}{12} \\ -\frac{1}{12} && -\frac{1}{12} && \frac{1}{4} - \lambda \\ \end{vmatrix} = 0 \\ \end{gather} $$
Определитель матрицы 3*3: $$ \begin{gather} (\frac{1}{4} - \lambda)^3 - 2 \cdot \frac{1}{12^3} - 3 \cdot \frac{1}{12^2} \cdot (\frac{1}{4} - \lambda) = 0 \\ \frac{(3-12\lambda)^3 - 2 - 3 \cdot (3 - 12 \lambda)}{12^3} = 0 \\ (3-12\lambda)^3 - 3 \cdot (3 - 12 \lambda) - 2 = 0 \\ \end{gather} $$
Пусть $t = 3 - 12 \lambda$ (т.е. $\lambda = \frac{3 - t}{12}$): $$ t^3 - 3 t - 2 = 0$$
Первый корень (перебором): $t_1 = -1 \to \lambda_1 = \frac{1}{3}$
$$ \begin{gather} (t + 1) (t^2 - t - 2) = 0 \\ t^2 - t - 2 = 0 \\ D = 1 + 8 = 9 \\ t_2 = \frac{1-3}{2} = -1 \to \lambda_2 = \frac{1}{3} \\ t_3 = \frac{1+3}{2} = 2 \to \lambda_3 = \frac{1}{12} \end{gather} $$
Собственные числа: $\frac{1}{12}$, $\frac{1}{3}$
Для каждого с.ч. найдем свой собственный вектор: $$v (A - \lambda I) = 0$$
- $\lambda_1 = \frac{1}{12}$:
$$ \begin{gather} \begin{cases} \frac{2}{12} x_1 -\frac{1}{12} x_2 -\frac{1}{12} x_3 = 0 \\ -\frac{1}{12} x_1 + \frac{2}{12} x_2 -\frac{1}{12} x_3 = 0 \\ -\frac{1}{12} x_1 -\frac{1}{12} x_2 + \frac{2}{12} x_3 = 0 \\ \end{cases} && \to v_1 = (1, 1, 1) \end{gather} $$
- $\lambda_2 = \frac{1}{3}$:
$$ \begin{gather} \begin{cases} -\frac{1}{12} x_1 -\frac{1}{12} x_2 -\frac{1}{12} x_3 = 0 \\ -\frac{1}{12} x_1 - \frac{1}{12} x_2 -\frac{1}{12} x_3 = 0 \\ -\frac{1}{12} x_1 -\frac{1}{12} x_2 - \frac{1}{12} x_3 = 0 \\ \end{cases} && \to v_2 = (-1, 1, 0), && v_3 = (-1, 0, 1) \end{gather} $$
Пары собственные числа-собственные значения: $$ \begin{gather} \lambda_1 = 1/12 && v_1 = (1, 1, 1) \\ \lambda_2 = 1/3 && v_2 = (-1, 1, 0) \\ \lambda_3 = 1/3 && v_3 = (-1, 0, 1) \end{gather} $$
- Проекция на главные компоненты
Главные компоненты - топ-K компонент, несущих бОльшую часть информацию (имеющих бОльшие собственные значения). С одной стороны, мы хотим снизить размерность данных, а с другой - сохранить больше информации.
Возьмем топ-$K$ компонент, которые кумулятивно объясняют ~90% вариации исходных данных. (На практике же, когда признаков очень много, число $K$ определяют "методом локтя"=elbow method)
Сначала отсортируем собственные вектора в соответствие с их собственными значениями:
$$\Lambda = (\lambda_2, \lambda_3, \lambda_1) = \big( \frac{1}{3}, \frac{1}{3}, \frac{1}{12} \big)$$
$$ V = (v_2, v_3, v_1) = \begin{pmatrix} -1 & -1 & 1 \\ 1 & 0 & 1 \\ 0 & 1 & 1 \\ \end{pmatrix} $$
И найдем долю объясненной дисперсии $d$, которую несет каждая компонента - по сути каждое собственное значение делим на сумму всех собственных значений.
$$ \begin{gather} d_2 = d_3 = \frac{1/3}{1/3+1/3+1/12} = \frac{4}{9} \\ d_1 = \frac{1/12}{1/3+1/3+1/12} = \frac{1}{9} \end{gather} $$
Чтобы понять, до какого числа компонент сокращать - должны посчитать долю дисперсии кумулятивно:
$$ \begin{gather} d_{PC_1} = 4/9 \\ d_{PC_2} = 8/9 \\ d_{PC_3} = 1 \end{gather} $$
Двумя главными компонентами можем объяснить 88.88 $\approx$ 90% дисперсии исходных данных. Поэтому третью компоненту скипаем, оставляем первые две. И далее делаем проекцию исходных данных на главные компоненты. Проецируя данные, мы переводим их из исходной системы координат (где переменные могут быть коррелированы) в новую, где каждая компонента независима от других.
Математически - ищем dot-product между исходными (центрированными) данными $Z$ и матрицей базисных векторов (топ-$K$ компонент) $V_K$: $$X' = Z \cdot V_K$$ $$ \begin{gather} X' = \begin{pmatrix} -0.25 & 0.75 & -0.25 \\ -0.25 & -0.25 & 0.75 \\ 0.75 & -0.25 & -0.25 \\ -0.25 & -0.25 & -0.25 \\ \end{pmatrix} \cdot \begin{pmatrix} -1 && -1 \\ 1 && 0 \\ 0 && 1 \\ \end{pmatrix} \\ X ' = \begin{pmatrix} 1 && 0 \\ 0 && 1 \\ -1 && -1 \\ 0 && 0 \end{pmatrix} \end{gather} $$
И код прилагается, но с.в. и трансформированная матрица будут различаться
# Initialize matrix
X = np.array([
[0, 1, 0],
[0, 0, 1],
[1, 0, 0],
[0, 0, 0]
])
# 1. Centralize matrix
means = X.mean(axis=0)
stds = np.sqrt(X.var(axis=0))
print(f'Means: {means.round(2)} \nStandard deviations: {stds.round(2)}')
Z = (X - means) # / stds
Z
Means: [0.25 0.25 0.25] Standard deviations: [0.43 0.43 0.43]
array([[-0.25, 0.75, -0.25],
[-0.25, -0.25, 0.75],
[ 0.75, -0.25, -0.25],
[-0.25, -0.25, -0.25]])
# 2. Compute covariance matrix
cov_matrix = np.cov(Z, rowvar=False)
cov_matrix
array([[ 0.25 , -0.08333333, -0.08333333],
[-0.08333333, 0.25 , -0.08333333],
[-0.08333333, -0.08333333, 0.25 ]])
# 3. Compute eigenvalues and eigenvectors
eigen_values, eigen_vectors = np.linalg.eig(cov_matrix)
print(f"Eigen-values: {eigen_values} \n\nEigen-vectors: \n{eigen_vectors.round(2)}")
Eigen-values: [0.33333333 0.08333333 0.33333333] Eigen-vectors: [[ 0.82 -0.58 -0.22] [-0.41 -0.58 -0.57] [-0.41 -0.58 0.79]]
# 4. Sort eigenvalues and eigenvectors in descending order of eigenvalues
order = np.argsort(eigen_values)[::-1]
eigen_vectors = eigen_vectors[:, order]
eigen_values = eigen_values[order]
# Find explained variance ratio to choose number of components
explained_variance_ratio_ = eigen_values / np.sum(eigen_values)
explained_variance_ratio_.cumsum().round(2)
array([0.44, 0.89, 1. ])
# 5. Project data into new space
n_components = 2
X_transformed = X.dot(eigen_vectors[:, :n_components])
X_transformed.round(4)
array([[-0.5709, -0.4082],
[ 0.791 , -0.4082],
[-0.2201, 0.8165],
[ 0. , 0. ]])
🧬 Task 2¶
Дано:
- Файлы с вариантами
chr21.pca.txt,chr22.pca.txt - Текстовый файл со списком сэплов
IBS.YRI.MEX.txt - Файл со списком сэплов и их принадлежностью к популяция
IBS.YRI.MEX.info.txt - Ноутбук
pop.model.l3.ipynbc PCA по 22ой хромосоме для популяций IBS, YRI, MEX.
def read_data(file, f_names, f_info):
with open(f_names, 'r') as f:
sample_names = [line.strip() for line in f]
with open(f_info, 'r') as f:
info_data = [line.strip().split('\t') for line in f]
info_dict = {row[0]: (row[1], row[3]) for row in info_data}
populations = [info_dict[name][1] for name in sample_names]
with open(file, 'r') as f:
genotypes = [
[
0 if allele == '0|0' else 2 if allele == '1|1' else 1
for allele in line.strip().split('\t')[3].split()
]
for line in f
]
genotype_matrix = np.array(genotypes).T
genotype_df = pd.DataFrame(genotype_matrix)
meta_df = pd.DataFrame({
'Pop': populations,
'ID': sample_names
})
result_df = pd.concat([genotype_df, meta_df], axis=1)
return result_df
Какой процент дисперсии описывает первая главная компонента для 22ой хромосомы?
gt22 = read_data(file=DATA_FOLDER / 'chr22.pca.txt',
f_names=DATA_FOLDER / 'IBS.YRI.MEX.txt',
f_info=DATA_FOLDER / 'IBS.YRI.MEX.info.txt'
)
print(gt22.shape)
gt22.head()
(278, 18271)
| 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | ... | 18261 | 18262 | 18263 | 18264 | 18265 | 18266 | 18267 | 18268 | Pop | ID | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 0 | 2 | 1 | 0 | 1 | 0 | 2 | 1 | 0 | 2 | ... | 0 | 1 | 0 | 0 | 0 | 1 | 0 | 0 | IBS | HG01503 |
| 1 | 1 | 1 | 0 | 1 | 1 | 1 | 1 | 1 | 0 | 2 | ... | 0 | 2 | 0 | 1 | 0 | 1 | 0 | 0 | IBS | HG01510 |
| 2 | 0 | 1 | 0 | 0 | 1 | 0 | 1 | 0 | 0 | 0 | ... | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | IBS | HG01515 |
| 3 | 0 | 2 | 0 | 0 | 2 | 0 | 2 | 0 | 1 | 1 | ... | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | IBS | HG01522 |
| 4 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | IBS | HG01527 |
5 rows × 18271 columns
X22 = gt22.drop(['Pop', 'ID'], axis=1)
pca22 = PCA(n_components=3)
X22_transformed = pca22.fit_transform(X22)
explained_variance_ratio_22 = pca22.explained_variance_ratio_
print(f'Первая компонента объясняет ~{100 * explained_variance_ratio_22[0]:.2f}% вариации данных')
Первая компонента объясняет ~15.64% вариации данных
plt.figure(figsize=(8, 6), dpi=80)
sns.barplot(x = [f'PC{x}' for x in range(1, 4)], y = explained_variance_ratio_22,
color='skyblue')
plt.xlabel('Principal components')
plt.ylabel('Explained variance ratio');
plt.figure(figsize=(8, 6), dpi=120)
sns.scatterplot(x=X22_transformed[:, 0], y=X22_transformed[:, 1], hue=gt22['Pop'].values)
plt.xlabel('PC1')
plt.ylabel('PC2')
plt.title("PCA on 22 chr genotype data");
Постройте PCA по 21ой хромосоме. Есть ли существенные отличия от визуализации по 22ой хромосоме?
gt21 = read_data(file=DATA_FOLDER / 'chr21.pca.txt',
f_names=DATA_FOLDER / 'IBS.YRI.MEX.txt',
f_info=DATA_FOLDER / 'IBS.YRI.MEX.info.txt'
)
print(gt21.shape)
gt21.head()
(278, 18860)
| 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | ... | 18850 | 18851 | 18852 | 18853 | 18854 | 18855 | 18856 | 18857 | Pop | ID | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 1 | 1 | 2 | 1 | 0 | 2 | 1 | 2 | 1 | 2 | ... | 1 | 1 | 1 | 1 | 2 | 1 | 0 | 2 | IBS | HG01503 |
| 1 | 0 | 1 | 2 | 1 | 1 | 1 | 1 | 1 | 2 | 1 | ... | 2 | 0 | 1 | 0 | 1 | 1 | 1 | 1 | IBS | HG01510 |
| 2 | 0 | 1 | 1 | 0 | 0 | 0 | 1 | 1 | 0 | 1 | ... | 2 | 0 | 1 | 0 | 2 | 2 | 0 | 2 | IBS | HG01515 |
| 3 | 1 | 1 | 1 | 1 | 0 | 1 | 1 | 1 | 1 | 1 | ... | 1 | 1 | 2 | 0 | 2 | 2 | 0 | 1 | IBS | HG01522 |
| 4 | 1 | 0 | 2 | 0 | 1 | 1 | 1 | 2 | 1 | 2 | ... | 2 | 0 | 1 | 0 | 1 | 1 | 1 | 1 | IBS | HG01527 |
5 rows × 18860 columns
X21 = gt21.drop(['Pop', 'ID'], axis=1)
pca21 = PCA(n_components=3)
X21_transformed = pca21.fit_transform(X21)
plt.figure(figsize=(8, 6), dpi=120)
sns.scatterplot(x=X21_transformed[:, 0], y=X21_transformed[:, 1], hue=gt21['Pop'].values)
plt.xlabel('PC1')
plt.ylabel('PC2')
plt.title("PCA on 21 chr genotype data");
Визуально отличий нет, YRI сгруппировались в кучу где-то отдельно (сильно отличаются от двух других популяций). MXL и IBS близки по первой компоненте, но разделяются второй. Также между ними есть небольшое пересечение
Спроектируйте наблюдения, построенных по 22ой хромосоме, на главные вектора, построенных по 21ой хромосоме, сравните с рисунком, где данные 22ой хромосомы спроектированы на вектора, полученных по 22ой хромосоме.
У 21 и 22 хромосом разное количество снипов. Не вижу другого выхода, кроме как обрезать по минимальному числу, иначе отображение не получится.
print(f'Features num: \nchr21 - {gt21.shape[1]}, chr22 - {gt22.shape[1]}')
Features num: chr21 - 18860, chr22 - 18271
num_features = min(X21.shape[1], X22.shape[1])
X21_filtered = X21.iloc[:, :num_features]
X22_filtered = X22.iloc[:, :num_features]
assert X21_filtered.shape == X22_filtered.shape
Переделываем PCA 21-й хромосомы, так как фичи изменились
pca21_new = PCA(n_components=3)
pca21_new.fit(X21_filtered)
PCA(n_components=3)In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.
PCA(n_components=3)
Проецируем наблюдения с 22-й хромосомы на вектора 21-й хромосомы
X22_transformed_from_21 = X22_filtered.dot(pca21_new.components_.T).to_numpy()
X22_transformed_from_21.shape
(278, 3)
Рисуем ✨картиночки✨
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12,6), dpi=120, sharex=True, sharey=True)
sns.scatterplot(x=X22_transformed[:, 0], y=X22_transformed[:, 1], hue=gt22['Pop'], ax=ax1)
sns.scatterplot(x=X22_transformed_from_21[:, 0], y=X22_transformed_from_21[:, 1], hue=gt22['Pop'], ax=ax2)
ax1.set_xlabel('PC1')
ax2.set_xlabel('PC1')
ax1.set_ylabel('PC2')
ax1.set_title('PCA on 22 chr')
ax2.set_title('PCA on 22 chr with 21 chr principal components');
Красота 💅
Объясняется тем, что PCA для 21-й хромосомы выделяет направления наибольшей дисперсии именно для снипов этой хромосомы. Главные компоненты 21-й хромосомы не отражают вариацию 22-й хромосомы, так как ее "ковариационная структура" выглядит иначе. Соответствия между важными позициями 21-й и 22-й хромосом, скорее всего нет (разве что на уровне случайности). Допустим, что в 21-й хромосоме важнее (более вариабельные) позиции {X}, а в 22-й - {Y}, и эти позиции не перекрываются. Главные компоненты для 21-й хромосомы будут переводить данные в пространство, где будет сохраняться вариабельность именно этих позиций. Если в 22-й хромосоме на этих позициях низкая вариабельность, или ее нет, то в новом пространстве ее и не станет. При этом для вариабельных позиций 22-й хромосомы при переходе в новое пространство исходня вариабельность затеряется, так как в матрице главных компонент 21-хромосомы этой структуры не было. В результате наблюдения 22-й хромосомы оказываются сжатыми в точку, так как по выбранным направлениям различий практически нет.
Объедините варианты 21ой и 22ой хромосом. Постройте PCA. Сравните результат с предыдущими.
X = pd.concat((X21, X22), axis=1)
X.columns = np.concatenate((X21.columns.values,
X22.columns.values + X21.columns.max()))
print(X.shape)
X.head()
(278, 37127)
| 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | ... | 37116 | 37117 | 37118 | 37119 | 37120 | 37121 | 37122 | 37123 | 37124 | 37125 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 1 | 1 | 2 | 1 | 0 | 2 | 1 | 2 | 1 | 2 | ... | 1 | 1 | 0 | 1 | 0 | 0 | 0 | 1 | 0 | 0 |
| 1 | 0 | 1 | 2 | 1 | 1 | 1 | 1 | 1 | 2 | 1 | ... | 2 | 0 | 0 | 2 | 0 | 1 | 0 | 1 | 0 | 0 |
| 2 | 0 | 1 | 1 | 0 | 0 | 0 | 1 | 1 | 0 | 1 | ... | 1 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 |
| 3 | 1 | 1 | 1 | 1 | 0 | 1 | 1 | 1 | 1 | 1 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 |
| 4 | 1 | 0 | 2 | 0 | 1 | 1 | 1 | 2 | 1 | 2 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 |
5 rows × 37127 columns
pca = PCA(n_components=3)
X_transformed = pca.fit_transform(X)
X_transformed.shape
(278, 3)
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(18, 6), dpi=120, sharex=True, sharey=True)
sns.scatterplot(x=X21_transformed[:, 0], y=X21_transformed[:, 1], hue=gt21['Pop'], ax=ax1)
sns.scatterplot(x=X22_transformed[:, 0], y=X22_transformed[:, 1], hue=gt22['Pop'], ax=ax2)
sns.scatterplot(x=X_transformed[:, 0], y=X_transformed[:, 1], hue=gt22['Pop'], ax=ax3)
for ax, name in zip((ax1, ax2, ax3), ('chr21', 'chr22', 'chr21 + chr22')):
ax.set_xlabel('PC1')
ax.set_ylabel('PC2')
ax.set_title(f'PCA on {name}')
plt.show()
Особо ничего не меняется. Только увеличилось "внутригрупповое" разнообразие внутри MXL - они стали более отдален друг от друга и от IBS (хотя некоторые наблюдения все равно пересекаются с этой популяцией)
Нарисуйте трехмерную проекцию на три главных компоненты для 22 хромосомы.
def plot_3d(x: np.ndarray, y: np.ndarray, z: np.ndarray,
populations: np.ndarray,
title: str = None, save_name: str = None
):
fig = px.scatter_3d(
x=x, y=y, z=z,
color=populations,
title=title,
opacity=0.8
)
fig.update_traces(marker=dict(size=4))
fig.update_layout(
scene=dict(
xaxis_title="PC1",
yaxis_title="PC2",
zaxis_title="PC3"
),
legend=dict(
title="Population",
x=1.05, y=1,
bgcolor="rgba(255, 255, 255, 0.8)"
),
margin=dict(l=0, r=0, b=0, t=30)
)
fig.show()
if save_name is not None:
fig.write_html(f"{save_name}.html")
plot_3d(x=X22_transformed[:, 0],
y=X22_transformed[:, 1],
z=X22_transformed[:, 2],
populations=gt22['Pop'].values,
title='PCA on chr22',
save_name='../imgs/HW1_task2_chr22'
)
Случайным образом уменьшите количество маркеров на 50 и на 80 процентов. Визуализируйте и проанализируйте.
Допустим это касается только 22-й хромосомы
# 1. Отбираем 50 и 20% случайных столбцов
X22_20 = X22.sample(frac=0.2, axis=1, random_state=SEED)
X22_50 = X22.sample(frac=0.5, axis=1, random_state=SEED)
X22.shape[1], X22_20.shape[1], X22_50.shape[1]
(18269, 3654, 9134)
# 2. Делаем pca
pca_22_20 = PCA(n_components=3)
pca_22_50 = PCA(n_components=3)
X22_20_transformed = pca_22_20.fit_transform(X22_20)
X22_50_transformed = pca_22_50.fit_transform(X22_50)
# 3. Визуализируем
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(18, 6), dpi=120, sharex=True, sharey=True)
sns.scatterplot(x=X22_transformed[:, 0], y=X22_transformed[:, 1], hue=gt22['Pop'], ax=ax1)
sns.scatterplot(x=X22_50_transformed[:, 0], y=X22_50_transformed[:, 1], hue=gt22['Pop'], ax=ax2)
sns.scatterplot(x=X22_20_transformed[:, 0], y=X22_20_transformed[:, 1], hue=gt22['Pop'], ax=ax3)
for ax, name in zip((ax1, ax2, ax3), ('chr22', 'chr22 (50% markers)', 'chr22 (20% markers)')):
ax.set_xlabel('PC1')
ax.set_ylabel('PC2')
ax.set_title(f'PCA on {name}')
plt.show()
Как при уменьшении до 50%, так и до 20% количества маркеров, заметно разделение популяций, хотя оно становится более "сжатым". Чем больше изначальных признаков мы удаляем, тем больше теряем информации о дисперсии - можем утратить важные направления, вдоль которых наблюдается вариация. Тем не менее, оказалось, что (особенно при 50%) это не так критично для описания различий между тремя популяциями, ведь они все равно различимы (хотя без раскраски это не было бы так очевидно, особенно для IBS и YRI при сокращении до 20% маркеров).
Структура данных сохраняется, вероятно, потому, что популяционные различия отражаются не отдельными маркерами а группами связанных маркеров. То есть многие маркеры несут избыточную информацию о популяционных различиях, поэтому удаление их части не критично, пока остаются маркеры, отражающие те же самые различия. Однако при сильном уменьшении количества маркеров вероятность того, что потеряются и такие группы, возрастает, что ведет к уменьшению "разрешающей способности" PCA
Уравняйте количество в каждой популяции. Изменится ли от этого PCA?
gt22.Pop.value_counts()
Pop YRI 108 IBS 106 MXL 64 Name: count, dtype: int64
Сделаем клиппинг по наименее распространенной популяции MXL. Скорее всего результат не изменится - не сильно выражена доминантность какой-то популяции (-ий).
pop_count = gt22.Pop.value_counts().min().item()
gt22_balanced = gt22.groupby('Pop').sample(pop_count, random_state=SEED)
X22_balanced = gt22_balanced.drop(['Pop', 'ID'], axis=1)
pca_22_balanced = PCA(n_components=3)
X22_balanced_transformed = pca_22_balanced.fit_transform(X22_balanced)
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12,6), dpi=120, sharex=True, sharey=True)
sns.scatterplot(x=X22_transformed[:, 0], y=X22_transformed[:, 1], hue=gt22['Pop'], ax=ax1)
sns.scatterplot(x=X22_balanced_transformed[:, 0], y=X22_balanced_transformed[:, 1], hue=gt22_balanced['Pop'], ax=ax2)
ax1.set_xlabel('PC1')
ax2.set_xlabel('PC1')
ax1.set_ylabel('PC2')
ax1.set_title('PCA on 22 chr')
ax2.set_title('PCA on 22 chr (balanced populations)');
Ничего существенно не изменилось, разве что цвета точек поменялись, что связано с изменившимся порядком в таблице.
Это указывает на то, что в изначальных данных преобладают межпопуляционные различия. То есть различия между популяциями выражены значительно сильнее, чем вариация внутри каждой популяции. Даже после удаления какого-то количества образцов из двух популяций (почти половина!) межгрупповая вариация сохраняется. Поскольку межгрупповая вариация доминирует над внутригрупповой, PCA ее захватывает и отображает в первых компонентах. Скорее всего, изменились только последние главные компоненты, которые отвечают за малые внутрипопуляционные различия.
🧬Task 3¶
Выберите 4 ваших любимых популяции из проекта 1000GP Phase 3.
Сделайте два текстовых файла: один со списком сэмплов, расположенных в столбик; второй состоит из двух столбцов
Sample,Population.Скачайте и запустите скрипт
pca.sh. Пропишите в нем путь до переменной VCF0 и путь до plink. Запустив его, вы получите текстовый файл, в котором данные отфильтрованы.Нарисуйте проекцию на две главные компоненты. Сопоставьте с тем, что известно про эти популяции из истории и географии.
BASE_URL = 'https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/'
VCF_URL = BASE_URL + 'ALL.chr22.phase3_shapeit2_mvncall_integrated_v5b.20130502.genotypes.vcf.gz'
ANNOTATION_URL = BASE_URL + 'integrated_call_samples_v3.20200731.ALL.ped'
VCF_PATH = DATA_FOLDER / '1KGP_phase3_22chr.vcf.gz'
SAMPLE_PATH = DATA_FOLDER / '1KGP_phase3_22chr_samples'
ANNOTATION_PATH = DATA_FOLDER / '1KGP_phase3_22chr_annotation'
- Скачиваем
.vcf22-й хромосомы и ее индекс
subprocess.run(f'wget {VCF_URL} -O {VCF_PATH} -q', shell=True)
subprocess.run(f'wget {VCF_URL}.tbi -O {VCF_PATH}.tbi -q', shell=True)
CompletedProcess(args='wget https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/ALL.chr22.phase3_shapeit2_mvncall_integrated_v5b.20130502.genotypes.vcf.gz.tbi -O ../data/HW1/1KGP_phase3_22chr.vcf.gz.tbi -q', returncode=0)
- Скачиваем список
sample_id, отбираем из них только те, которые соответствуют выбранным популяциям- GBR - British in England and Scotland
- CHB - Han Chinese in Beijing, China
- LWK - Luhya in Webuye, Kenya
- JPT - Japanese in Tokyo, Japan
Решила усложнить себе жизнь
Еще, чтобы воспользоваться функцией чтения файлов, я делала файл не sample-популяция, а sample-gender-sample-population (я не поняла, что в 3 колонке должно быть, но она и не используется)
samples_script ="""
#!/bin/bash
DATA_FOLDER=$HOME/PopGen_Models/data/HW1
ANNOTATION_URL="https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/integrated_call_samples_v3.20200731.ALL.ped"
ANNOTATION_PATH=$DATA_FOLDER/1KGP_phase3_22chr_annotation
SAMPLE_PATH=$DATA_FOLDER/1KGP_phase3_22chr_samples
# Get all samples
bcftools query -l ${DATA_FOLDER}/1KGP_phase3_22chr.vcf.gz > all_samples
# Install annotation
wget ${ANNOTATION_URL} -q -O tmp.csv
# Query populations -> 'SAMPLE-POPULATION' file
awk '$7 == "GBR" || $7 == "CHB" || $7 == "LWK" || $7 == "JPT" {print $2 "\t" $5 "\t" $2 "\t" $7}' tmp.csv > annotation_samples
# Intersect with .vcf samples
join -1 1 -2 1 -t $'\t' <(sort all_samples) <(sort annotation_samples) > ${ANNOTATION_PATH}
echo "Annotation saved in ${ANNOTATION_PATH}"
# 'SAMPLE' file
awk '{print $1}' ${ANNOTATION_PATH} > ${SAMPLE_PATH}
echo "Samples saved in ${SAMPLE_PATH}"
rm tmp.csv annotation_samples all_samples
"""
SAMPLE_EXEC_PATH = SCRIPT_FOLDER / 'get_samples.sh'
with open(SAMPLE_EXEC_PATH, 'w') as f:
print(samples_script, file=f)
subprocess.run(f'chmod +x {SAMPLE_EXEC_PATH}', shell=True)
subprocess.run(f'bash {SAMPLE_EXEC_PATH}', shell=True)
Annotation saved in /home/ksu_marshmallow/PopGen_Models/data/HW1/1KGP_phase3_22chr_annotation Samples saved in /home/ksu_marshmallow/PopGen_Models/data/HW1/1KGP_phase3_22chr_samples
CompletedProcess(args='bash ../scripts/get_samples.sh', returncode=0)
- Запускаем скрипт
pca.sh. Я в нем только изменила пути, чтобы ко мне в папку с данными все скачивалось + изменила название выходного файла, чтобы не перезаписывался файл из второй таски.
PCA_EXEC_PATH = SCRIPT_FOLDER / 'pca.sh'
subprocess.run(f'chmod +x {PCA_EXEC_PATH}', shell=True)
subprocess.run(f'{PCA_EXEC_PATH} 22 {SAMPLE_PATH}', shell=True)
PLINK v1.90b6.21 64-bit (19 Oct 2020) www.cog-genomics.org/plink/1.9/ (C) 2005-2020 Shaun Purcell, Christopher Chang GNU General Public License v3 Logging to /home/ksu_marshmallow/PopGen_Models/data/HW1/test2.log. Options in effect: --allow-extra-chr --double-id --indep-pairwise 50 10 0.1 --out /home/ksu_marshmallow/PopGen_Models/data/HW1/test2 --set-missing-var-ids @:# --vcf /home/ksu_marshmallow/PopGen_Models/data/HW1/test.vcf.gz 15374 MB RAM detected; reserving 7687 MB for main workspace. --vcf: /home/ksu_marshmallow/PopGen_Models/data/HW1/test2-temporary.bed + /home/ksu_marshmallow/PopGen_Models/data/HW1/test2-temporary.bim + /home/ksu_marshmallow/PopGen_Models/data/HW1/test2-temporary.fam written. 1054487 variants loaded from .bim file. 1054487 missing IDs set. 397 people (0 males, 0 females, 397 ambiguous) loaded from .fam. Ambiguous sex IDs written to /home/ksu_marshmallow/PopGen_Models/data/HW1/test2.nosex . Using 1 thread (no multithreaded calculations invoked). Before main variant filters, 397 founders and 0 nonfounders present. Calculating allele frequencies... 10111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697989 done. 1054487 variants and 397 people pass filters and QC. Note: No phenotypes present. Pruned 770911 variants from chromosome 22, leaving 283576. Pruning complete. 770911 of 1054487 variants removed. Marker lists written to /home/ksu_marshmallow/PopGen_Models/data/HW1/test2.prune.in and /home/ksu_marshmallow/PopGen_Models/data/HW1/test2.prune.out . PLINK v1.90b6.21 64-bit (19 Oct 2020) www.cog-genomics.org/plink/1.9/ (C) 2005-2020 Shaun Purcell, Christopher Chang GNU General Public License v3 Logging to /home/ksu_marshmallow/PopGen_Models/data/HW1/test3.log. Options in effect: --allow-extra-chr --double-id --extract /home/ksu_marshmallow/PopGen_Models/data/HW1/test2.prune.in --maf 0.1 --make-bed --out /home/ksu_marshmallow/PopGen_Models/data/HW1/test3 --set-missing-var-ids @:# --vcf /home/ksu_marshmallow/PopGen_Models/data/HW1/test.vcf.gz 15374 MB RAM detected; reserving 7687 MB for main workspace. --vcf: /home/ksu_marshmallow/PopGen_Models/data/HW1/test3-temporary.bed + /home/ksu_marshmallow/PopGen_Models/data/HW1/test3-temporary.bim + /home/ksu_marshmallow/PopGen_Models/data/HW1/test3-temporary.fam written. 1054487 variants loaded from .bim file. 1054487 missing IDs set. 397 people (0 males, 0 females, 397 ambiguous) loaded from .fam. Ambiguous sex IDs written to /home/ksu_marshmallow/PopGen_Models/data/HW1/test3.nosex . --extract: 283576 variants remaining. Using 1 thread (no multithreaded calculations invoked). Before main variant filters, 397 founders and 0 nonfounders present. Calculating allele frequencies... 10111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697989 done. 266392 variants removed due to minor allele threshold(s) (--maf/--max-maf/--mac/--max-mac). 17184 variants and 397 people pass filters and QC. Note: No phenotypes present. --make-bed to /home/ksu_marshmallow/PopGen_Models/data/HW1/test3.bed + /home/ksu_marshmallow/PopGen_Models/data/HW1/test3.bim + /home/ksu_marshmallow/PopGen_Models/data/HW1/test3.fam ... 101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899done.
CompletedProcess(args='../scripts/pca.sh 22 ../data/HW1/1KGP_phase3_22chr_samples', returncode=0)
- Делаем PCA
gt22_task3 = read_data(file=DATA_FOLDER / 'task3_chr22.pca.txt',
f_names=SAMPLE_PATH,
f_info=ANNOTATION_PATH
)
print(gt22_task3.shape)
gt22_task3.head()
(397, 17186)
| 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | ... | 17176 | 17177 | 17178 | 17179 | 17180 | 17181 | 17182 | 17183 | Pop | ID | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 1 | 1 | 0 | 1 | 2 | 1 | 1 | 0 | 0 | 0 | ... | 1 | 1 | 0 | 1 | 0 | 0 | 1 | 0 | GBR | HG00096 |
| 1 | 0 | 1 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 1 | ... | 1 | 1 | 1 | 1 | 1 | 0 | 1 | 0 | GBR | HG00097 |
| 2 | 0 | 1 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 1 | ... | 0 | 1 | 0 | 0 | 0 | 1 | 0 | 0 | GBR | HG00099 |
| 3 | 0 | 2 | 1 | 0 | 1 | 0 | 2 | 0 | 1 | 2 | ... | 1 | 0 | 0 | 0 | 0 | 2 | 0 | 0 | GBR | HG00100 |
| 4 | 0 | 2 | 1 | 0 | 1 | 0 | 2 | 1 | 1 | 2 | ... | 1 | 1 | 0 | 1 | 0 | 2 | 0 | 0 | GBR | HG00101 |
5 rows × 17186 columns
# Распределение популяций
gt22_task3.Pop.value_counts()
Pop JPT 104 CHB 103 LWK 99 GBR 91 Name: count, dtype: int64
X22_task3 = gt22_task3.drop(['Pop', 'ID'], axis=1)
pca22_task3 = PCA(n_components=3)
X22_task3_transformed = pca22_task3.fit_transform(X22_task3)
plt.figure(figsize=(8, 6), dpi=120)
sns.scatterplot(x=X22_task3_transformed[:, 0], y=X22_task3_transformed[:, 1], hue=gt22_task3['Pop'].values)
plt.xlabel('PC1')
plt.ylabel('PC2')
plt.title("PCA on 22 chr genotype data");
Давным-давно, в самой древней Африке, жили люди. Это было время, когда все люди на Земле были как одна большая семья. Они жили рядом друг с другом и знали только один дом - африканские равнины.
Однажды эта семья решила отправиться в большое путешествие. Кто-то остался дома, в теплых и солнечных краях (это были предки наших кенийцев, LWK), а другие захотели увидеть мир. И вот они начали свой долгий путь...
Некоторые из них пошли на северо-запад, через горы и пустыни. Дорога была нелегкой, но они нашли новые земли - холодные и суровые. Эти смельчаки стали предками европейцев, а потом и британцев (GBR). В этих новых местах им пришлось учиться жить в холоде, одеваться теплее и находить новую еду. Их тела и гены начали меняться, чтобы приспособиться к новому (дивному) миру.
А другие отправились совсем в другую сторону - на восток. Они шли через Индию, через горы и реки, пока не добрались до великой Азии. Там они тоже разделились: одни остались в Китае (CHB), а другие через какое-то время переплыли море и обосновались на островах Японии (JPT). Хотя они теперь жили далеко друг от друга, их предки когда-то были одной семьей, поэтому их гены все еще очень похожи.
Те, кто остался в Африке (LWK), продолжали жить так же, как их предки. Они знали каждый уголок своей родной земли, приспособились к ее жаре и дождям. Их гены оставались особенными, потому что они никогда не покидали свой дом.
Спустя тысячи лет, когда ученые начали изучать гены людей со всего мира, они заметили удивительную вещь. Гены тех, кто живет в Европе (британцы), совсем не похожи на гены тех, кто живет в Африке (кенийцы). А гены китайцев и японцев, хотя они и разные, все еще напоминают друг друга - ведь их предки когда-то были одной семьей.
И вот, когда мудрые ученые построили свою волшебную карту PCA, она показала три большие семьи: африканцы (LWK) остались в своем древнем доме, европейцы (GBR) ушли далеко на север, а азиаты (CHB и JPT) отправились на восток. Первые две главные дороги на карте (PC1 и PC2) рассказали о большом разделении между континентами. А китайцы с японцами, хотя и живут отдельно уже давно, все еще стоят близко — их история расхождения спрятана где-то на более тонких тропинках, может быть, на третьей компоненте.
Так генетика помогла нам прочитать древнюю историю человечества, записанную в наших клетках 🌍✨
plot_3d(x=X22_task3_transformed[:, 0],
y=X22_task3_transformed[:, 1],
z=X22_task3_transformed[:, 2],
populations=gt22_task3['Pop'].values,
title='PCA on chr22',
save_name='../imgs/HW1_task3_chr22'
)
Действительно! Третья компонента PCA уже может разделять японцев и китайцев. Хоть тоже не идеально и они пересекаются, но хотя бы не накладываются друг на друга, как только с первыми двумя векторами.